#include "DissectBubblesCorrector.h"
#include <queue>
#include <stack>
#include <map>
#include <iostream>
#include <cmath>
using namespace std;

int DissectBubblesCorrector::indel_penalty = 1;
int DissectBubblesCorrector::replace_penalty = 1;

SequenceGraph* SimpleMergeCorrector::correctGraph()
{
	SequenceVertex *vertex;
	const vector<SequenceEdge> *pred,*succ;
	m_seqGraph->resetCursor();
	while (m_seqGraph->hasNext())
	{
		vertex = m_seqGraph->next();
		succ = &(vertex->getSuccessors());
		while (succ->size()==1 && (*succ)[0].end->numPredecessors() == 1)
		{
			if (m_seqGraph->mergeVertexIntoPred(vertex,(*succ)[0].end)) {
				++m_numCorrected;
				succ = &(vertex->getSuccessors());
			}	else {
				break;
			}
		} // while (succ->size()==1 && (*succ)[0]->numPredecessors() == 1)
		pred = &(vertex->getTwinVertex()->getSuccessors());
		while (pred->size() == 1 && (*pred)[0].end->numPredecessors() == 1)
		{
			if (m_seqGraph->mergeVertexIntoSucc((*pred)[0].end->getTwinVertex(),vertex)) {
				++m_numCorrected;
				pred = &(vertex->getTwinVertex()->getSuccessors());
			}	else {
				break;
			}
		} // while (pred->size() == 1 && (*pred)[0]->numPredecessors() == 1)
	} // while (m_seqGraph->hasNext())
	return m_seqGraph;
}

void DissectBubblesCorrector::TraceRecord::record(GID_T seqVID, int treeID, TraceNode* node) {
	iterator_type findIter = m_map.find(seqVID);
	if (findIter == m_map.end()) {
		vector<vector<TraceNode*> > val;
		for (int i = 0; i < treeCount; ++i) {
			val.push_back(vector<TraceNode*>());
		}
		findIter = (m_map.insert(make_pair(seqVID, val))).first;
	} // if (findIter == m_map.end())
	(findIter->second)[treeID].push_back(node);
}


vector<pair<int, DissectBubblesCorrector::TraceNode*> > DissectBubblesCorrector::TraceRecord::otherPossiblePaths(const TraceNode& node, int treeID, int absMismatch, int pruneBound)
{
	vector<pair<int, TraceNode*> > ret;
	iterator_type findIter = m_map.find(node.seqVertexID);
	if (findIter == m_map.end()) {
		return ret;
	}
	for (int i = 0; i < treeCount; ++i) {
		if (i == treeID) { continue; }
		for (vector<TraceNode*>::iterator iter = (findIter->second)[i].begin(); iter != (findIter->second)[i].end();) {
			if ((*iter)->dist < pruneBound) {
				iter = (findIter->second)[i].erase(iter); // prune
			} else {
				if (absMismatch > (node.dist - (*iter)->dist)) {
					ret.push_back(make_pair(i, *iter));
				}
				++iter;
			}
		}
	}
	return ret;
}

void DissectBubblesCorrector::TraceRecord::removeRecord(int treeID)
{
	for (iterator_type it = m_map.begin(); it != m_map.end(); ++it) {
		(it->second)[treeID].clear();
	}
}

DissectBubblesCorrector::TraceRecord::~TraceRecord()
{
	for (iterator_type iter = m_map.begin(); iter != m_map.end(); ++iter) {
		iter->second.clear();
	}
}

vector<DissectBubblesCorrector::TraceNode*> DissectBubblesCorrector::DistNodeMap::vicinities( int dist, int absMismatch, int treeID )
{
	VLOG(2) << "Enter vicinities()";
	vector<TraceNode*> ret;
	int nearDist = dist - absMismatch;
	int farDist = dist + absMismatch;
	iterator_type lowerIter = m_map.lower_bound(nearDist);
	iterator_type upperIter = m_map.upper_bound(farDist);

	for (iterator_type iter = lowerIter; iter != upperIter; ++ iter) {
		for (int i = 0; i < iter->second.size(); ++i) {
			if (i != treeID) {
				ret.insert(ret.end(), (iter->second)[i].begin(), (iter->second)[i].end());
			}
		} // for (int i = 0; i < iter->second.size(); ++i)
	} // for (iterator_type iter = m_map.begin(); iter != m_map.end();)

	for (iterator_type iter = m_map.begin(); iter != lowerIter; ++iter) {
		iter->second.clear();
	}
	m_map.erase(m_map.begin(), lowerIter);

	VLOG(2) << "Exit vicinities()";
	return ret;
}

void DissectBubblesCorrector::DistNodeMap::insert( int treeID, TraceNode* node )
{
	int dist = node->dist;
	iterator_type iter = m_map.find(dist);
	if (iter == m_map.end()) {
		vector<vector<TraceNode*> > val;
		for (int i = 0; i < m_treeCnt; ++i) {
			val.push_back(vector<TraceNode*>());
		}
		iter = (m_map.insert(make_pair(dist, val))).first;
	} // if (iter == m_map.end())
	(iter->second)[treeID].push_back(node);
}

DissectBubblesCorrector::DistNodeMap::~DistNodeMap()
{
	for (iterator_type iter = m_map.begin(); iter != m_map.end(); ++iter)
	{
		iter->second.clear();
	}
}

DissectBubblesCorrector::DissectBubblesCorrector(SequenceGraph* seqGraph,
	int searchBound, int mismatchUpperBound, int mismatchLowerBound, double relativeMismatch):
		AbstractGraphCorrector(seqGraph), m_searchBound(searchBound),
		m_mismatchUpperBound(mismatchUpperBound), m_mismatchLowerBound(mismatchLowerBound), m_relativeMismatch(relativeMismatch),
		m_simpleMerger(seqGraph)
{
	int arrSize = (seqGraph->maxVertexID() >> 3) + 1;
	m_visitRecord = new unsigned char[arrSize];
	m_traceRecord = new unsigned char[arrSize];
	memset(m_visitRecord, 0, arrSize);
	memset(m_traceRecord, 0, arrSize);
	distMatrix = new int[m_searchBound << 1];
}

DissectBubblesCorrector::~DissectBubblesCorrector()
{
	delete[] m_visitRecord;
	delete[] distMatrix;
}


SequenceGraph* DissectBubblesCorrector::correctGraph()
{
	cout <<  "Start dissect bubbles..." <<endl;
	LOG(INFO) <<  "Start dissect bubbles...";
	SequenceVertex* root = NULL;
	m_seqGraph->resetCursor();
	while (root = nextRoot()) {
		VLOG(2) << "next root: " <<root->getVertexID();
		dissectBubble(root);
	} // while (root = nextRoot())
	m_simpleMerger.correctGraph();
	cout <<"Merged Errors: " <<m_simpleMerger.numCorrected() <<endl;
	cout <<  "Finish dissect bubbles." <<endl;
	LOG(INFO) <<  "Finish dissect bubbles.";
	return m_seqGraph;
}

SequenceVertex* DissectBubblesCorrector::nextRoot()
{
	VLOG(2) << "Enter nextRoot()";
	while (m_seqGraph->hasNext()) {
		SequenceVertex* ret = m_seqGraph->next();
		VLOG(2) <<"considering " <<ret->getVertexID();
		if (!isVisited(ret->getVertexID()) && m_seqGraph->numOfPredecessors(ret) == 0) {
			VLOG(2) <<"Exit nextRoot() with ret";
			return ret;
		} else 	if (!isVisited(ret->getVertexID()) && m_seqGraph->numOfSuccessors(ret) == 0) {
			VLOG(2) <<"Exit nextRoot() with ret";
			return ret->getTwinVertex();
		}
	} // if (m_seqGraph->hasNext())
	VLOG(2) <<"Exit nextRoot() with NULL";
	return NULL;
}

void DissectBubblesCorrector::dissectBubble(SequenceVertex* root)
{
	VLOG(2) << "Enter dissectBubble";
	if (isVisited(root->getVertexID())) {
		return;
	}
	queue<SequenceVertex*> q;
	setVisited(root->getVertexID());
	q.push(root);
	while (!q.empty()) {
		SequenceVertex* curV = q.front();
		q.pop();
		const vector<SequenceEdge> succs(m_seqGraph->getSuccessors(curV)); // this is a snapshot
		const vector<SequenceEdge>& succs_active = m_seqGraph->getSuccessors(curV);
		for (int i = 0; i < succs.size(); ++i) {
			if (!isVisited(succs[i].end->getVertexID())) {
				setVisited(succs[i].end->getVertexID());
				q.push(succs[i].end);
			} else if (!isTraced(succs[i].end->getVertexID()) ) {
				if (!contain(succs[i], succs_active)) { continue; }
				setTraced(succs[i].end->getVertexID());
				tryResolve(succs[i].end);
			}
		} // for (int i = 0; i < succs.size(); ++i)
	} // while (!q.empty())
	VLOG(2) << "Exit dissectBubble";
}

//bool DissectBubblesCorrector::isTerminated(const vector<FibHeap>& queues, int tid)
//{
//	VLOG(2) << "Enter isTerminated()";
//	if (queues[tid].isEmpty()) { return false; }
//	if (((TraceNode*) queues[tid].peekMin())->dist == INT_MAX) {
//		VLOG(2) << "Exit isTerminated()";
//		return true;
//	} // if (((TraceNode*) queues[tid].peekMin())->dist == INT_MAX)
//	VLOG(2) << "Exit isTerminated()";
//	return false;
//}

void DissectBubblesCorrector::tryResolve(SequenceVertex* end) {
	VLOG(2) << "Enter tryResolve()";
	// see which sequence node is traced on what nodes on what trees
	vector<SequenceEdge> preds(m_seqGraph->getPredecessors(end));
 	int treeCnt = preds.size();
	vector<bool> detached(treeCnt, false); // if the tree is already detached from end
	for (int i = 0; i < treeCnt - 1; ++i) {
		for (int j = i + 1; j < treeCnt; ++j) {
			if (!detached[i] && !detached[j] && preds[i].multiplicity != preds[j].multiplicity) { //multiplicity same is impossible to detach
				set<GID_T> candidates;
				if (generateCandidates(preds[i].end, preds[j].end, end, candidates)) {
					tryResolve(preds, i, j, end, detached, candidates);
				}
			} // if (!detached[i] && !detached[j] && multiplicity[i] != multiplicity[j])
		} // for (int i = 0; i < treeCnt - 1; ++i) {for (int j = i + 1; j < treeCnt; ++j)
	} // for (int i = 0; i < treeCnt - 1; ++i)
	VLOG(2) << "Exit tryResolve()";
}

void DissectBubblesCorrector::tryResolve(const vector<SequenceEdge>& allpreds, int one, int another, SequenceVertex* end, vector<bool>& detached, const set<GID_T>& candidates)
{
	VLOG(1) << "Enter tryResolve() level2";
	const int limit = 1000000;
	int cnt = 2;
	vector<SequenceVertex*> preds;
	preds.push_back(allpreds[one].end);
	preds.push_back(allpreds[another].end);
	vector<FibHeap> queues(2, FibHeap(TraceNode::cmp_func));
	vector<TraceNode> roots(2, TraceNode(0, NULL));
	TraceRecord traceRecord(2);
	DistNodeMap distMap(2);

	for (int i = 0; i < 2; ++i) {
		roots[i].seqVertexID = preds[i]->getVertexID();
		//roots[i].parent = NULL;
		roots[i].dist = 0;
		if ( preds[i]->size() - KmerVertex::getLength() + 1 <= m_searchBound) {
			// must be sure elements of roots then never change,
			// since roots[i]'s address could change due to expansion
			queues[i].insert(&(roots[i]));
		}
		traceRecord.record(roots[i].seqVertexID, i, &(roots[i]));
		distMap.insert(i, &(roots[i]));
	} // for (int i = 0; i < treeCnt; ++i)
	while (true) {
		int curTreeId = chooseLowestTree(queues); // id of lowest tree
		if ( curTreeId == -1 || ((TraceNode*)queues[curTreeId].peekMin())->dist >= m_searchBound) { // ensure path will not exceed m_searchBound
			VLOG(1) << ((curTreeId == -1) ? "queue empty" : "search bound reached");
			break;
		}
		TraceNode* curNode = (TraceNode*)queues[curTreeId].extractMin();
		//cerr <<"dequeue dist: " << curNode->dist <<endl;
		SequenceVertex* curV = m_seqGraph->getVertexByID(curNode->seqVertexID);
		if ( curNode->dist && (!curV->numPredecessors()
			|| !shouldGrow(curNode, curTreeId, distMap)) ) { continue; }

		const vector<SequenceEdge>& curPreds = m_seqGraph->getPredecessors(curV);
		curNode->children = vector<TraceNode>(curPreds.size(), TraceNode(0, curNode));
		for (int i = 0; i < curPreds.size(); ++i) {
			curNode->children[i].seqVertexID = curPreds[i].end->getVertexID();
			//curNode->children[i].parent = curNode;
			curNode->children[i].dist = curNode->dist + curV->size() - KmerVertex::getLength() + 1;
			if ( !candidates.count(curNode->children[i].seqVertexID) || cnt >= limit) { continue; }
			//cerr <<"new item enqueue dist: " <<curNode->children[i].dist <<endl;
			traceRecord.record(curNode->children[i].seqVertexID, curTreeId, &(curNode->children[i]) );
			distMap.insert(curTreeId, &(curNode->children[i]));
			const vector<pair<int, TraceNode *> >& otherPossiblePaths =
				traceRecord.otherPossiblePaths(curNode->children[i], curTreeId, m_mismatchUpperBound, curNode->dist + 1 - m_mismatchUpperBound);
			// the last parameter should be curNode->dist + 1 - absmismatch
			// curNode->dist + 1 is the minimal length of path could be encounter later
			if (otherPossiblePaths.empty()) { // if not empty it should not grow
				if (curNode->children[i].dist + curPreds[i].end->size() - KmerVertex::getLength() + 1 <= m_searchBound) { // ensure path will not exceed m_searchBound
					queues[curTreeId].insert(&(curNode->children[i])); // risky code since you have to guarantee that children never expand
					++cnt;
				}
			} // if (otherPossiblePaths.empty())
			for (int j = 0; j < otherPossiblePaths.size(); ++j) {
				int otherTreeId = otherPossiblePaths[j].first;
				//if (isTerminated(queues, otherTreeId)) { continue; }
				//TraceNode& collide = curNode->children[i];
				TraceNode* np1 = &(curNode->children[i]);
				TraceNode* np2 = otherPossiblePaths[j].second;
								
				// use multiplicity to prune impossible detach
				SequenceVertex *begin = m_seqGraph->getVertexByID(np1->seqVertexID);
				SequenceVertex *succ1 = np1->parent ? m_seqGraph->getVertexByID(np1->parent->seqVertexID) : end;
				SequenceVertex *succ2 = np2->parent ? m_seqGraph->getVertexByID(np2->parent->seqVertexID) : end;
				SequenceVertex *pred1 = m_seqGraph->getVertexByID( roots[curTreeId].seqVertexID );
				SequenceVertex *pred2 = m_seqGraph->getVertexByID( roots[otherTreeId].seqVertexID );
				//SequenceVertex *end is pass-in parameter
				int beginPlicity1 = m_seqGraph->getMultiplicity(begin, succ1);
				int beginPlicity2 = m_seqGraph->getMultiplicity(begin, succ2);
				int endPlicity1 = m_seqGraph->getMultiplicity(pred1, end);
				int endPlicity2 = m_seqGraph->getMultiplicity(pred2, end);
				if ( !((beginPlicity1 > beginPlicity2 && endPlicity1 > endPlicity2) ||
					(beginPlicity1 < beginPlicity2 && endPlicity1 < endPlicity2)) ) {
						continue;
				}

				vector<SequenceVertex*> path1; // current tree path
				vector<SequenceVertex*> path2; // rival tree path
				while ((np1 = np1->parent) != NULL)
				{
					path1.push_back(m_seqGraph->getVertexByID(np1->seqVertexID));
				} // while (np1 != NULL)
				while ((np2 = np2->parent) != NULL)
				{
					path2.push_back(m_seqGraph->getVertexByID(np2->seqVertexID));
				} // while (np1 != NULL)
				// recover 2 paths and compare
				recoverPath(path1, p1content);
				recoverPath(path2, p2content);
				int aThreshold = max(ceil(p2content.size() * m_relativeMismatch), ceil(p1content.size() * m_relativeMismatch));
				if (areSimilar(p1content, p2content, max(aThreshold, m_mismatchLowerBound) )) {
					// detach untrusted one
					if (beginPlicity1 > beginPlicity2	&& endPlicity1 > endPlicity2) {
							m_seqGraph->removeSuccessorAndTwinPredecessor(begin, succ2);
							++m_numCorrected;
							if (begin != pred2 || end != succ2) {
								m_seqGraph->removeSuccessorAndTwinPredecessor(pred2, end);
							}
							// record detach
							detached[otherTreeId == 0 ? one : another] = true;
							cleanTrees(roots);
							VLOG(1) << "Exit tryResolve() level2 with resolve";
							return;
							// terminate the tree
							//queues[otherTreeId].clear();
							//roots[otherTreeId].dist = INT_MAX;
							//queues[otherTreeId].insert(&roots[otherTreeId]);
							//traceRecord.removeRecord(otherTreeId);
					} // if (m_seqGraph->getMultiplicity(begin, succ1) > m_seqGraph->getMultiplicity(begin, succ2) && m_seqGraph->getMultiplicity(pred1, end) > m_seqGraph->getMultiplicity(pred2, end))
					else if (beginPlicity1 < beginPlicity2	&& endPlicity1 < endPlicity2) {
							m_seqGraph->removeSuccessorAndTwinPredecessor(begin, succ1);
							++m_numCorrected;
							if (begin != pred1 || end != succ1) {
								m_seqGraph->removeSuccessorAndTwinPredecessor(pred1, end);
							}
							// record detach
							detached[curTreeId == 0 ? one : another] = true;
							cleanTrees(roots);
							VLOG(1) << "Exit tryResolve() level2 with resolve";
							return;
							//// terminate the tree
							//queues[curTreeId].clear();
							//roots[curTreeId].dist = INT_MAX;
							//queues[curTreeId].insert(&roots[curTreeId]);
							//traceRecord.removeRecord(curTreeId);
							//break;
					} // if (m_seqGraph->getMultiplicity(begin, succ2) > m_seqGraph->getMultiplicity(begin, succ1) && m_seqGraph->getMultiplicity(pred2, end) > m_seqGraph->getMultiplicity(pred1, end))
				} // if (areSimilar(p1content, p2content))
			} // for (int j = 0; j < otherPossiblePaths.size(); ++j)
			//if (isTerminated(queues, curTreeId)) {
			//	break;
			//}
		} // for (int i = 0; i < curPreds.size(); ++i)
	} // while (true)	
	cleanTrees(roots);
	VLOG(1) << "Exit tryResolve() level2 without resolve";
}

int DissectBubblesCorrector::generateCandidates(SequenceVertex* pred1, SequenceVertex* pred2, SequenceVertex* end, set<GID_T>& candidates)
{
	VLOG(1) << "Enter generateCandidate()";
	const int limit = 1000000;
	int cnt = 0;
	vector<CandidateNode> comEnds(2, CandidateNode(end->getVertexID(), 0));
	comEnds[1].branch = 1;
	vector<SequenceVertex*> preds;
	preds.push_back(pred1);
	preds.push_back(pred2);
	map<GID_T, CandidateNode> visited;
	FibHeap minHeap(CandidateNode::cmp_func);
	for (int i = 0; i < 2; ++i) {
		CandidateNode candi(preds[i]->getVertexID(), i);
		candi.parents.push_back(&(comEnds[i]) );
		minHeap.insert(&(visited.insert(make_pair(candi.seqVertexID, candi)).first->second));
		++cnt;
	} // for (int i = 0; i < 2; ++i)

	int searchBound = m_searchBound;
	CandidateNode* common = NULL;
	while (!minHeap.isEmpty() && ((CandidateNode*)minHeap.peekMin())->dist < searchBound) {
		CandidateNode* curNode = static_cast<CandidateNode*>(minHeap.extractMin());
		SequenceVertex* curVertex = m_seqGraph->getVertexByID(curNode->seqVertexID);
		int nextDist = curNode->dist + curVertex->size() - KmerVertex::getLength() + 1;
		if (nextDist > searchBound) { continue; }
		int curBranch = curNode->branch;
		const vector<SequenceEdge>& curPreds = curVertex->getPredecessors();
		for (int i = 0; i < curPreds.size(); ++i) {
			map<GID_T, CandidateNode>::iterator it = visited.find(curPreds[i].end->getVertexID());
			if (it != visited.end() ) {
				if (!common && it->second.branch != curBranch) {
					common = &(it->second);
					searchBound = min(common->dist + m_mismatchUpperBound, m_searchBound);
				}
				(it->second).parents.push_back(curNode);
			} else {
				CandidateNode candi(curPreds[i].end->getVertexID(), curBranch);
				candi.dist = nextDist;
				candi.parents.push_back(curNode);
				if (++cnt <= limit) {
					minHeap.insert(&(	visited.insert(make_pair(candi.seqVertexID, candi)).first->second));
				}
			}
		} // for (int i = 0; i < curPreds.size(); ++i)
	} // while (!minHeap.isEmpty() && ((CandidateNode*)minHeap.peekMin())->dist < searchBound)
	
	if (common == NULL) {
		VLOG(1) << "Exit generateCandidate() with 0";
		return 0;
	}
	vector<int> group[2];
	for (int i = 0; i < common->parents.size(); ++i) {
		CandidateNode* candi = common->parents[i];
		group[candi->branch].push_back(m_seqGraph->getMultiplicity(common->seqVertexID, candi->seqVertexID));
	} // for (int i = 0; common->parents.size(); ++i)
	bool flag = false; // indicate if a detach is ever possible
	for (int i = 0; i < group[0].size() && !flag; ++i) {
		for (int j = 0; j < group[1].size(); ++j) {
			if (group[0][i] != group[1][j]) { flag = true; break; }
		}
	} // for (int i = 0; i < group[0].size(); ++i)
	if (flag) {
		queue<CandidateNode*> q;
		q.push(common);
		while (!q.empty()) {
			CandidateNode* candi = q.front();
			candidates.insert(candi->seqVertexID);
			q.pop();
			for (int i = 0; i < candi->parents.size(); ++i) {
				if (!candidates.count((candi->parents)[i]->seqVertexID) )
				q.push((candi->parents)[i]);
			}
		} // while (!q.empty()) 
	}
	VLOG(1) << "Exit generateCandidate() with " << candidates.size();
	return candidates.size();
}

bool DissectBubblesCorrector::contain(const SequenceEdge& ele, const vector<SequenceEdge>& collection) const
{
	for (int i = 0; i < collection.size(); ++i) {
		if (ele == collection[i]) {
			return true;
		}
	} // for (int i = 0; i < collection.size(); ++i)
	return false;
}

int DissectBubblesCorrector::chooseLowestTree( const vector<FibHeap>& queues )
{
	VLOG(2) << "Enter chooseLowestTree()";
	int ret = -1;
	int minHeight = INT_MAX;
	for (int i = 0; i < queues.size(); ++i) {
		if (!queues[i].isEmpty()) {
			int d = ((TraceNode*) queues[i].peekMin())->dist;
			if (d < minHeight) {
				ret = i;
				minHeight = d;
			} // if (d < minHeight)
		} // if (!queues[i].isEmpty())
	} // for (int i = 0; i < queues.size(); ++i)
	VLOG(2) <<"Exit chooseLowestTree()";
	return ret;
}

int DissectBubblesCorrector::chooseLowestBranch( const vector<FibHeap>& queues )
{
	VLOG(2) << "Enter chooseLowestBranch()";
	int ret = -1;
	int minHeight = INT_MAX;
	for (int i = 0; i < queues.size(); ++i) {
		if (!queues[i].isEmpty()) {
			int d = ((CandidateNode*) queues[i].peekMin())->dist;
			if (d < minHeight) {
				ret = i;
				minHeight = d;
			} // if (d < minHeight)
		} // if (!queues[i].isEmpty())
	} // for (int i = 0; i < queues.size(); ++i)
	VLOG(2) <<"Exit chooseLowestBranch()";
	return ret;
}

void DissectBubblesCorrector::recoverPath(const vector<SequenceVertex*>& pathVertices, BitDeque& pathContent)
{
	VLOG(2) << "Enter recoverPath()";
	pathContent.clear();
	int pos = 0;
	for (int i = 0; i < pathVertices.size(); ++i) {
		int numKmer = pathVertices[i]->size() - KmerVertex::getLength() + 1;
		if (pos + numKmer < KmerVertex::getLength()) {
			pos += numKmer;
			continue;
		} // if (pos + pathVertices[i]->m_numOfKmer < KmerVertex::getLength())
		int start = KmerVertex::getLength() - 1 - pos;
		pos +=  numKmer;
		if (start < 0) {
			pathContent.push_back_range(pathVertices[i]->m_compressedSequence, 0, pathVertices[i]->size() - (KmerVertex::getLength() - 1));
		} else {
			pathContent.push_back_range(pathVertices[i]->m_compressedSequence, start, pathVertices[i]->size() - (KmerVertex::getLength() - 1));
		}
	} // for (int i = 0; i < pathVertices.size(); ++i)
	VLOG(2) <<"Exit recoverPath()";
}

bool DissectBubblesCorrector::areSimilar(const BitDeque& str1, const BitDeque& str2, int threshold)
{
	VLOG(2) << "Enter areSimilar()";
	int len1 = str1.size();
	int len2 = str2.size();
	// below is not needed here since it has already checked before
	//if (abs(len1 - len2) > m_mismatchUpperBound) {
	//	return false;
	//}
	if (len1 == 0) {
		return (len2 * indel_penalty <= threshold);
	}
	if (len2 == 0) {
		return (len1 * indel_penalty <= threshold);
	}

	memset(distMatrix, 0, (m_searchBound << 1) * sizeof(int));
	for (int i = 0; i < m_searchBound; ++i) {
		distMatrix[i] = (i + 1) * indel_penalty;
	}
	int* row0 = distMatrix;
	int* row1 = distMatrix + m_searchBound;
	for (int i = 0; i < len1; ++i) {
		int rowMin = INT_MAX;
		for (int j = 0; j < len2; ++j) {
			if (j == 0) {
				row1[j] = min(i + ((str1[i] == str2[j]) ? 0 : replace_penalty), min((i + 1) * indel_penalty + indel_penalty, row0[j] + indel_penalty));
				rowMin = min(row1[j], rowMin);
			} else {
				row1[j] = min(row0[j - 1] + ((str1[i] == str2[j]) ? 0 : replace_penalty), min(row1[j -1] + indel_penalty, row0[j] + indel_penalty));
				rowMin = min(row1[j], rowMin);
			}
		} // {for (int j = 0; j < len2; ++j)
		if (rowMin > threshold) {
			VLOG(2) << "Exit areSimilar() with false";
			return false;
		}
		swap(row0, row1);
	} // for (int i = 0; i < len1; ++i)
	if (row0[len2 - 1] <= threshold) {
		VLOG(2) << "Exit areSimilar() with true";
		return true;
	}
	VLOG(2) << "Exit areSimilar() with false";
	return false;
}

int DissectBubblesCorrector::getArrIdx(char ch) {
	switch (ch) {
		case 'A': return 0;
		case 'C': return 1;
		case 'G': return 2;
		case 'T': return 3;
		default: LOG(FATAL) << "Unknown base";
	}
}

bool DissectBubblesCorrector::shouldGrow( TraceNode* curNode, int curTreeId, DistNodeMap& distMap )
{
	VLOG(2) << "Enter shouldGrow()";
	const int limit = 10;
	const vector<TraceNode*>& vicinities = distMap.vicinities(curNode->dist, m_mismatchUpperBound, curTreeId);
	for (int i = 0; i < vicinities.size(); ++i) {
		vector<SequenceVertex*> path1, path2;
		TraceNode* np1 = curNode;
		TraceNode* np2 = vicinities[i];
		int estimatedLen1 = np1->dist;
		int estimatedLen2 = np2->dist;
		while ((np1 = np1->parent) != NULL) {
			path1.push_back(m_seqGraph->getVertexByID(np1->seqVertexID));
		}
		while ((np2 = np2->parent) != NULL) {
			path2.push_back(m_seqGraph->getVertexByID(np2->seqVertexID));
		}

		if (path1.size() >= limit || path2.size() >= limit) {
			return false;
		}
		recoverPathCompleteHead(path1, p1content);
		recoverPathCompleteHead(path2, p2content);

		///*cerr <<"p1: ";
		//for (int ii = 0; ii < path1.size(); ++ii) {
		//	cerr <<path1[ii]->getVertexID() <<' ';
		//}
		//cerr <<endl;
		//cerr <<"p2: ";
		//for (int ii = 0; ii < path2.size(); ++ii) {
		//	cerr <<path2[ii]->getVertexID() <<' ';
		//}
		//cerr <<endl;
		//cerr <<p1content <<endl;
		//cerr <<p2content <<endl;*/

		if (estimatedLen2 >= estimatedLen1) { // compare loosely
			int mismat = ceil(estimatedLen2 * m_relativeMismatch);
			if (areSimilar(p1content, p2content, min(mismat, m_mismatchUpperBound << 1))) {
				VLOG(2) << "Exit shouldGrow with true.";
				return true;
			}
		} else { // compare tighly
			int mismat = ceil(estimatedLen1 * m_relativeMismatch);
			if (areSimilar(p1content, p2content, min(mismat, m_mismatchUpperBound))) {
				VLOG(2) << "Exit shouldGrow with true.";
				return true;
			}
		}
	} // for (int i = 0; i < vicinities.size(); ++i)
	VLOG(2) << "Exit shouldGrow with false.";
	return false;
}

void DissectBubblesCorrector::recoverPathCompleteHead(const vector<SequenceVertex*>& pathVertices, BitDeque& pathContent)
{
	VLOG(2) << "Enter recoverPathCompleteHead()";
	pathContent.clear();
	for (int i = 0; i < pathVertices.size(); ++i) {
		pathContent.push_back_range(pathVertices[i]->m_compressedSequence, 0, pathVertices[i]->size() - (KmerVertex::getLength() - 1));
	} // for (int i = 0; i < pathVertices.size(); ++i)
	VLOG(2) <<"Exit recoverPathCompleteHead()";
}

void DissectBubblesCorrector::cleanTrees( vector<TraceNode>& roots )
{
	stack<TraceNode*> stk;
	for (int i = 0; i < roots.size(); ++i) {
		stk.push(&(roots[i]));
	}
	bool traversed = false;
	while (!stk.empty()) {
		TraceNode* cur = stk.top();
		if (traversed) {
			cur->children.clear();
			stk.pop();
			traversed = false;
			if (stk.empty() || stk.top() == cur->parent) {
				traversed = true;
			}
		} else if (!cur->children.size()) {
			stk.pop();
			if (stk.empty() || stk.top() == cur->parent) {
				traversed = true;
			}
		} else {
			for (int i = 0; i < cur->children.size(); ++i) {
				stk.push(&(cur->children[i]));
			}
		}
	} // while (!stk.empty())
}